%% ===========================[description] ============================
%---------------- implements the function-----------------
% With multiple inputs
% Draw amplitude and frequency       
%---------------- prefiles -----------------
%        quasi0StiffSuspensionCertain.slx
%        quasi0StiffSuspension.slx
%        n1nDpolyRoots.m
%        GetStatistics.m/GetStat.m
% ------------------ enter -------------------
% Cylinder air pressure
% Piston distance from rear wall
% Cylinder node distance
% Effective area of cylinder
%% ========================[Clear workspace]============================
clear
close all;
%% =========================[parameter area] ===============================
d0Certain=32.5e-3;                     %piston distance from rear wall m
AcCertain=pi*0.1^2;                    %cylinder piston area m^2
L0Certain=0.467;                       %cylinder node distance m
kCertain=3;                            %cylinder air pressure atm
%--------------------------------------------------------------------------
mCertain=800;
m=800;                                 % mass on spring kg
m1=100;                                % tyre mass
Pa=101*10^3;
A_level=16;                            %-|
B_level=64;                            %-|-road
C_level=256;                           %-|
Level=B_level;
v=60;                                  %Vehicle speed
C=2000;                                %damping
qA=0.1;                                %Excitation amplitude
fM=0.1:0.1:30;
lengthoffre=length(fM);
%-------------------------------------------------------------------------
excSwitch=-1;                          %Select sine excitation
D=length(1);                           %number of random variables
order=3;
coeOfVar=[0.1,0,0,0];
vMat=[d0Certain,AcCertain,L0Certain,kCertain];
Collocation=n1nDpolyRoots(order,D); 
groupNumber=length(Collocation(:,1));
PCmean=zeros(lengthoffre,3);
PCvar=zeros(lengthoffre,3);
dqrms=zeros(lengthoffre,1);
%% ==========================[Computing]================================
for i=1:lengthoffre
    f=fM(i);
    for groNum=0:groupNumber  
        if groNum==0
            [d0,Ac,L0,k]=iniRanVar(d0Certain,AcCertain,L0Certain,kCertain,...
                         coeOfVar,Collocation,groupNumber,1);
            sim('quasi0StiffSuspension.slx')
            outputData=zeros(3,groupNumber);
        else
            [d0,Ac,L0,k]=iniRanVar(d0Certain,AcCertain,L0Certain,kCertain,...
                         coeOfVar,Collocation,groupNumber,groNum);
            sim('quasi0StiffSuspension.slx')
            outputData(:,groNum)=[rms(a.data),rms(fd.data),rms(Fd.data)];
        end
    end
    for j=1:3
        [PCmean(i,j),PCvar(i,j)]=GetStat ...
          (Collocation,outputData(j,:)',order);
    end
    dqrms(i)=rms(dq.data);
    
    %----------------------------------------------------------------------
    proPers=floor(i/lengthoffre*100);                                    %
    clc                                                                  %----Progress bar
    disp(['PCE[','#'*ones(1,proPers),' '*...                             %
        ones(1,100-proPers) ,']',num2str(proPers),'%'])                  %
    %----------------------------------------------------------------------
end
Hout=PCmean+2*sqrt(PCvar);
Lout=PCmean-2*sqrt(PCvar);
aqpH=Hout./dqrms;
aqpL=Lout./dqrms;
%% =============================[Drawing]============================
close all
yl={'$$\left| \frac{\ddot{z}}{\dot{q}}\right|\rm{(s^{-1})}$$','$$\left|\frac{f_d}{\dot{q}}\right|\rm{(s)}$$','$$\left| \frac{F_d}{G\dot{q}}\right|\rm{(s/m)}$$'};
for p=1:3
    figure(p)
    loglog(fM,aqpH(:,p),fM,aqpL(:,p),'LineWidth',1.5)
    set(gca,'ygrid','on','gridlinestyle','--','Gridalpha',0.8,'MinorGridalpha',0.5)
    grid on;
    xlabel('\it{Frequence}\rm{(Hz)}')
    ylabel(yl(p),'Interpreter','latex','FontSize',14)
    legend('Maximum Boundary','Minimum Boundary')
end




